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Abstract 

Hadronic matrix elements of operators relevant to nucleon decay in grand unified theories are 

(N 

' calculated numerically using lattice QCD. In this context, the domain-wall fermion formulation, 

(N ' 

^> • combined with non-perturbative renormalization, is used for the first time. These techniques bring 

o , 

\ reduction of a large fraction of the systematic error from the finite lattice spacing. Our main effort 

is devoted to a calculation performed in the quenched approximation, where the direct calculation 
of the nucleon to pseudoscalar matrix elements, as well as the indirect estimate of them from the 
JL . nucleon to vacuum matrix elements, are performed. First results, using two flavors of dynamical 

(— i \ domain-wall quarks for the nucleon to vacuum matrix elements are also presented to address the 



X 



systematic error of quenching, which appears to be small compared to the other errors. Our 
results suggest that the representative value for the low energy constants from the nucleon to 
vacuum matrix elements are given as \a\ ~ \(3\ ~ 0.01 GeV 3 . For a more reliable estimate of the 
physical low energy matrix elements, it is better to use the relevant form factors calculated in the 
direct method. The direct method tends to give smaller value of the form factors, compared to 
the indirect one, thus enhancing the proton life-time; indeed for the ir° final state the difference 
between the two methods is quite appreciable. 

PACS numbers: 12.38.Gc, 11.30.Rd, 12.10.Dm 
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I. INTRODUCTION 



In the standard model baryon number is conserved up to a very good approximation. 
While it is broken (only) by the electro-weak anomaly, the size of breaking is extremely 
small so that baryon number violation would be extremely difficult to observe in any high 
energy experiment. The conservation of baryon number happens to be exact in the classical 
level in the standard model. However, it is not a consequence of the underlying fundamental 
symmetries, e. g. gauge and Lorentz symmetries. In general, baryon number is not protected 
once we extend the gauge groups of the standard model to larger symmetry groups. The 
existence of the baryon number violating interaction makes protons unstable. This particular 
phenomena, proton decay, or in general, nucleon decay is one of the most decisive signals 
that prove that nature has larger fundamental symmetry than that of the standard model. 
Grand unified theories (GUTs) 0>Q|, with or without supersymmetry, possess such a feature. 

Although proton decay has been searched for in the deep mine experiments for more than 
a decade, it has not been observed j^j]. The lower bound of the proton lifetime from the 
experiments has provided an important stringent constraint on the GUT models. In fact, 
the simplest minimal SUSY GUT with SU(5) gauge group is almost surely excluded 3| |s8 1 . 

At low energies, the processes that allow nucleon decay may be represented in terms 
of a low energy effective Hamiltonian made up of Standard Model fields. This effective 
Hamiltonian will be dominated by the operators of lowest dimension. The requirement that 
the low energy effective theory has the same symmetry as the standard model is strong 
enough to constrain the form of the possible operators in the effective Hamiltonian 0, |(| [J 
The operators consist of three quarks and one lepton fields as (q c q)(J C q), which lead to a 
decay of a nucleon to a pseudoscalar and an anti-lepton. As the low energy hadrons are 
involved in the decay, this is a highly non-perturbative process, and so lattice QCD provides 
an ideal tool to analyze it. 

Since difficult to calculate the hadronic matrix element (PS\(q c q)q\N) of the three-quark 
operator with initial nucleon (N) and final pseudoscalar (PS) states, various QCD model 
approaches were used to calculate nucleon to vacuum matrix element (0\(q c q)q\N) , which 
in turn gives (PS\(q c q)q\N) with the help of chiral perturbation theory. In 80's many 

n nn n n n nn 

calculations of the matrix elements [8|, |9|, HO, [11|, [12J, [13|, HJ, [15|| were reported. These 
calculations varied by a factor of O(10) between smallest to largest estimate (see and 
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Table IVl]) . As partial width of the decay is proportional to the matrix elements squared, 
different model calculations lead to a factor of 0(100) difference in proton lifetime. 

nnn 

161 1171 1181 ] using lattice QCD were performed in late 80 's. These 



The first calculations 



calculations also took the approach of calculating the nucleon to vacuum matrix elements, 
using chiral perturbation theory to infer the final result. The nucleon to vacuum matrix 
elements obtained in these studies disagreed as well, with a factor of about 5 difference 
between the tallest Q to largest Q e 8tim ate. 

There are various sources of systematic errors involved in such calculations, one of which 
is the large range extrapolation of the matrix elements in the quark mass. As the matrix 
elements are used as the low energy constants of the chiral Lagrangian, they must be obtained 
in the chiral limit. These original lattice studies necessarily used unphysically heavy values 
of the up and down quark masses. In fact, the values used were typically of the order of the 
mass of the physical strange quark. 

Another source of error comes from measuring the lattice scale. Since the low energy 
constants measured have mass dimension three, their error - when quoted in physical units 
- receives a contribution equal to three times the error on this scale. There are various ways 
of extracting the lattice scale which need not agree at finite lattice spacing, or when working 
in the quenched approximation. Moreover, some conventional ways to set the lattice scale 
such as using the p meson or the nucleon mass as an input again, require a large extrapolation 
in the mass to reach the physical point. It is important to note that all the calculations 
performed up to now used unimproved Wilson fermions; a formulation known to have large 
scaling violations. It is not surprising that a moderate difference in the lattice cut off leads 
to a large difference in the low energy constants. 

Operator renormalization is another source of systematic error. The operators which 
mediate the nucleon decay must be renormalized in some renormalization scheme. This 
must be the same renormalization scheme that was used for the calculation of the Wilson 
coefficients for the effective Hamiltonian, and is usually MS scheme. The common way 
to perform this calculation has been to use lattice perturbation theory, which has poor 
convergence property, leading to a large systematic error. 

There is one important systematics error, which neither the model or lattice calculations 
in the 80 's address: even if one gets the correct value of the low energy constants, the matrix 
element for the decay could have an appreciable error. The reason is that the pion resulting 
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from the decay has, in the center of mass frame, a momentum around half of the nucleon 
mass, where the leading order chiral expansion may not be a good approximation. 

The, more recent, JLQCD work 3] was better at addressing many of these systematics 
compared to the old results: they employed smaller quark masses, and used an improved 
lattice perturbation theory. Most importantly, the first reliable direct measurement for the 
nucleon to pseudoscalar matrix elements were performed. However, the calculation was still 
performed in the quenched approximation, and made use of the unimproved Wilson action. 
Later, a joint collaboration including some of the original JLQCD members addressed the 
issue of the scaling violations in the restricted case of the nucleon to vacuum matrix elements 
Q. It appeared the value that JLQCD g ot at a finite lattice S paci„ g was .a rg et by ahnost 
factor 2 than that in the continuum limit, showing the particular difficulty in taking the 
continuum limit with the unimproved Wilson fermions. 

In this paper we calculate the nucleon decay matrix elements in a approach similar to 
JLQCD, but, with three key differences: we use domain-wall fermions, which preserve chiral 
symmetry to a very high accuracy (and so are expected to have much reduced scaling 
violations), we use a non-perturbative renormalization scheme, and we also investigate the 
effects of moving away from the quenched approximation. As mentioned previously, all 
the works which have been done so far for nucleon decay used the quenched approximation 
wherein the quark loops in the propagation of the gluons are neglected. Although it has been 
shown that light flavored hadron spectrum is reproduced by the quenched approximation 
within 10% difference of experiment [2]], the systematic error due to this approximation is 
process dependent and, in general, uncontrollable. We present the first study of quenching 
error by calculating the nucleon to vacuum matrix elements using unquenched u, d quarks 
and comparing to the quenched result. As a result of these differences an appreciable 
reduction of the systematic error is expected. 

This paper is organized as follows: in Sect. |H] we summarize the general properties of the 
baryon number violating operator and the nucleon decay matrix elements. Sect. IIHI discusses 
operator property in view of the renormalization and gives the detail of how to apply the 
non-perturbative renormalization scheme to the nucleon decay operator. The calculation of 
nucleon decay matrix elements with the quenched approximation is shown in Sec. IIVI The 
estimate of the nucleon to vacuum matrix elements with unquenching u, d quark is given in 
Sec. El Sec. EH is devoted to the conclusions. Some results obtained in earlier stages of this 
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work have been reported in Refs. 




We use the Euclidean lattice formulation, so the metric and gamma matrices should be 
taken as Euclidean. Dimensional quantities are often written in the lattice units (a). We 
avoid writing "a" explicitly in most cases. 



II. GENERAL DESCRIPTION 



In this section we describe the operators and matrix elements to calculate and their 
theoretical background. Lattice calculation setup is overviewed, which in the later sections 
will be discussed in detail. 



A. Properties of operators and matrix elements 

The general type of the nucleon decay operator that appears in the low energy effective 
Hamiltonian is constrained by the symmetries of the standard model, SU(3) C x SU(2)l x 
U(1)y j^llf]]- The quark part must contain three quark fields, as the SU(3) C singlet made 
of 3 x 3 x 3 is the lowest dimensional baryon number violating operator which is invariant 
under SU (3) c . Another fermion field (lepton or antilepton) is necessary to make the operator 
Lorentz invariant. The SU(2)l x U(1)y symmetry rules out the possibility of having the 
antilepton field in the operator, ajid restricts the flavor/chirality combination to four types 

these operators read 



antilepton field in the operator, and restricts the flavor Jd 
In the notation of Weinberg Q]; Abbott and Wise p] 



0$L = (DlUl) R (g^,f d ) L e^ (1) 

0$L = {CA%(Pll d ) R e^ (2) 

<52L = {CA%^A)Le ijk e^e^ (3) 

0$L = (DlUiUU^h)^, (4) 

where I is the generic lepton field, q is the left handed quark field, U and D denote 
up and down type right handed quarks. a,b,c,d are generation numbers, k are 
SU(3) color labels, a,fl,j,8 are SU(2) indices. The inner product (x,y)R/L is defined as 
( X ,y)n/L = x T CPu/Ly, where C is the charge conjugation matrix and Pr/l is the right/left 
handed projection matrix. The vector and tensor Dirac matrices have been eliminated in 
the expression by Fierz rearrangement. These dimension six operators [3] are the lowest 



dimensional operators that appear in the low energy effective Hamiltonian written in terms 
of Standard Model particles. Higher dimensional operators are suppressed by inverse powers 
of the heavy mass scale characteristic of the fundamental high energy theory (eg. Mx for 
GUTs). As is evident from the form of the operator, it breaks baryon number (B), but pre- 
serves baryon minus lepton number (B — L), leading to a decay of nucleon to a pseudoscalar 
and an antilepton. 

The general form of the three quark part of the operator, which transforms as a spinor, 

is 

OS = (ud) r s r = e ilk (u iT CP v d j )P r s\ (5) 

where r(r') can be either R or L. u, d, and s are quark fields not necessarily labeling the real 
u, d, s flavors. From the form of the operator it is evident that the proton or neutron cannot 
decay to a final state that has strangeness S < 0. Thus, processes such as p —> K° + l + 
are disallowed. The trivial constraint that the mass of the pseudoscalar (PS) is below that 
of the nucleon (N), nips < m^, requires the final state to be one of 7r°, tt^, K°, K + or 77. 
Hence the real physical flavor of a quark in Eq. © is one of the three lightest, u, d, s. 

We calculate the hadronic matrix element with the nucleon and the allowed pseudoscalar 
states, 

(PS;p\O rr '\N;k,s), (6) 

where p and k are spatial momenta of pseudoscalar and nucleon respectively, and s — ±1/2 
is the spin of the nucleon. Parity transformation yields a relation between different chirality 
matrix elements: 

(PS;p\0 R/LL \N;k,s) =^(PS;-p\0 L/RR \N;-k,s). (7) 

For the range of accuracy that is expected from our calculation, it is sufficient to assume 
isospin symmetry, which further reduces the number of independent matrix elements. Fol- 
lowing Ref. Q| here we list 14 matrix elements and their "isospin partners" obtained by 
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exchanging u and d: 



(ir \(ud) R / L u L \p) = 


(7i \(du) R/L d L \n), 


(8) 


(n + \(ud) R/L d L \p) = 


-(ir~\(du) R/L u L \n), 


(9) 


(K°\(us)b/tUt\v) = 
\i\ \\u,o ) R/L U 'L\1 J / 


-(K + \(ds) D /rdr \n) 


do) 


(K + \(us) R/L d L \p) = 


-(K°\(ds) R/L u L \n), 


(11) 


(K + \(ud) R/L s L \p) = 


-(K \(du) R/L s L \n), 


(12) 


(K + \(ds) R/L u L \p) = 


-(K°\(us) R/L d L \n), 


(13) 


(v\(ud) R/L u L \p) = 


-(r)\{du) R/L d L \n). 


(14) 



Our interpolating field for each hadron state is summarized in appendix [XJ The negative 
signs on the rhs of Eqs. (J%j)- (fi4*j) appear from the transformation of the interpolating fields 
J n o — > —J^o and J p — > —J n by interchanging u and d. There is a relation between final tt° 
and final tt^ matrix elements in the isospin limit 

(K + \(ud) R/L d L \p) = V2(7r°\(ud) R/L u L \p). (15) 

We call the Ihs of Eqs. (jHJ), (fTUjl - ffT^ the principal matrix elements. In the quenched sim- 
ulation, we are going to calculate these 12 principal matrix elements. All the other pos- 
sible matrix elements are obtained from the principal matrix elements by Eqs. (|%|)- ()14j) . 
Eq. ©and Eq. |TK|) . Note, however, flavor SU(3) breaking effect of rj is not treated in this 
paper [7l|. A nucleon to pseudoscalar decay is characterized by its initial three momentum 
k, spin s = ±1/2 and final momentum p. By Lorentz covariance the matrix element is 

n 

required to have the form [19], 

(PS;p\0 R / LL \N; k, s) = P L [W R/LL (q 2 ) - iiWf' L L (q 2 )}u N (k, s), (16) 

where is the nucleon spinor. The form factors Wq and W q are functions of the square 
of the momentum transfer = k^ — p^. For the physical decay this is the momentum 
of the lepton. However, on the lattice we work with unphysical values of the masses and 
momentum transfer, then extrapolate to the physical point. In this case, Wq and W q are the 
functions of m^, mps and q 2 . For the range of masses and momenta used in our simulation 
the two terms in the square brackets seem to be of the same order. This fact implies that 
the second term is negligible in the physical amplitude due to the on-shell condition on the 
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lepton (— q 2 



mr 







. Hence, Wq is called relevant and W„ is irrelevant form factor. As 



mr 



0, Wo(0) is the final target of our calculation. The parity condition Eq. (J2J) 



implies the relation between form factors with different chirality as 



W R/L 



R( J2\ 



W^ RL {q% 



(17) 



for x = 0, q. Note that parity holds only after the statistical average over lattice gauge field 
configurations. On a single gauge configuration the rhs and Ihs of Eq. ()17j) generally differ, 
and we can take average of them to get better statistics, which is done in our analysis. 
Once the form factor is calculated the partial width of the p — ^ PS + 1 decay is obtained 

as 



r(p —> ps + 7) 



m r 



32tt 2 



m PS 



m r 



^CW^p^PS) 



(18) 



where m p (mps) denotes the proton (pseudoscalar) mass, and we have set the lepton mass to 
zero. C l are the Wilson coefficients of the dimension-six operator in the low energy effective 
Hamiltonian, 

£ 6 = Y,C%q)(ql)Y = ~J2° l W C O m ] 1 ■ (19) 

i i 

The index i distinguishes the type (flavor and chirality) of the three quark operator O qqq , 
which is one of those in the matrix elements Eqs. (f%|)- ()14)) . as well as the type of lepton. 
While Cj and Wq depends on the renormalization scheme and scale, the dependence cancels 
out in the product. In this work we use a renormalization scale \i = 2 GeV and MS scheme 
with the naive dimensional regularization (NDR) for Wq. With Wq calculated in this work, 
and a knowledge of Cj, which depend on the particular GUT model, proton lifetime can be 
estimated through Eq. |TR|) . 



B. Usage of chiral perturbation theory 

In this work we discuss the calculation of the form factors for the principal matrix elements 
by means of numerical lattice simulations. To obtain the results at the physical kinematics in 
terms of quark masses and pseudoscalar momentum, it is necessary to know the dependence 
of the form factors on these parameters, especially for the u, d quark masses as these are not 
attainable on the lattice with the present algorithmic and computational resources. Chiral 
perturbation theory (x?T) gives such an information. The tree level result for the form 



8 



factors of the nucleoli decay are available jisj . The "direct method", which uses three point 
functions and various two-point functions, estimates the matrix elements for the kinematics 
of the particular lattice simulation. Then they are extrapolated/interpolated to the physical 
kinematics with the help of x?T. 

The chiral Lagrangian of the nucleon decay j^j] involves only two additional low energy 
parameters to the ordinary three-flavor baryon chiral Lagrangian at leading order. Measuring 
these parameters on the lattice and combining with the other parameters of the baryon chiral 
Lagrangian, all the matrix elements of nucleon decay can, in principle, be calculated. For 
the proton to 7r° decay as an example, the relevant form factors read 

W RL (p -> vr ) = a(l + D + F)/V2f, (20) 
W LL (p -> vr ) = (3(1 + D + F)/V2f, (21) 

where / is the tree level pion decay constant with a normalization such that the experimental 
value is f n ~ 131 MeV. D and F are the couplings of baryons to the axial current, where 
the sum of them gives the nucleon axial charge: D + F = g^. a and (3 are specific to 
the nucleon decay which can be calculated at leading order through the proton to vacuum 
matrix element of the operators, 

aP L u p = (0|Q2S,b>, PPl»p = <0|aL L >>, (22) 

aP R u p = -(0\O£*\p), (3P R u p = -(0\0™\P), (23) 



where Eq. (|23|) is obtained from Eq. (j22|) by parity transformation. We fix the phase defini- 
tion such that a and (3 are real and a < 0. As we will later describe, we observe a + (3 ~ 0, 
which is expected because of the relation, 

(a + (3)u p = -(0\e^ k (u Tl Cd^ l5 u k \p), (24) 

which vanishes in the non-relativistic limit and is known to be quite small even at small 
quark maS8 e 8 Q. 

Reduction formulae similar to Eqs. fl2()L (1211) are available for all the principal matrix 
elements These are summarized in Appendix iBl This way of calculating the matrix 

elements is referred to as the "indirect method" . It has to be noted that the indirect method 
can have sizable systematic error, which is difficult to estimate reliably: 
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1. In the £77(3)/ baryon chiral Lagrangian, there are four parameters which control the 
flavor SU (3) breaking effects. Two parameters are used to match the baryon masses, 
and enter the nucleon decay matrix elements through baryon masses. The other two 
contribute for the nucleon to pseudo-scalar matrix element we wish to extract (and 
calculate in the direct method), but do not contribute to the nucleon to vacuum 
matrix elements required in the indirect method. As we have no means to estimate 
their values, we set these parameters, named bi (i — 1,2) in Refs. [tqI . 25 1, as bi — 



as is standard for such calculations. However, they can naturally be 0(1). Setting 
\b{\ = 1 for a test, the contribution of the breaking term to the relevant form factors 
is estimated as 5-30% for N — * K decays. The worst case is (r]\(ud) R UL\p) , where the 



contribution is as large as 200% 
should be negligible. 



73j . For N — > it the effect of these parameters 



2. Even if flavor violation effects were small, i. e., effects of bi were negligible, there could 
be appreciable systematic error from the usage of the lowest order x?T. This is due 
to the large energy of the pseudoscalar, which is at least a half of the nucleon mass: 
Eps > m^jl in the CM frame, while lowest order %PT is exact only at the zero 
energy (soft pion) limit. Of course, the direct method also relies upon leading order 
chiral perturbation theory at an energy scale around m^jl. This is clearly a source 
of systematic error. However, since the extrapolation required is over a much shorter 
distance, it may be expected the systematic error is smaller than that for the indirect 
method. 

The indirect method requires the calculation of only a few two-point functions on the 
lattice. The direct method is superior to the indirect method since for the former there is 
no need to assume any parameter in the chiral Lagrangian to be in any particular range 
of values. Also, as mentioned, the former has less reliance on The coefficients are 

determined by fit to the lattice results obtained for each decay process independently. The 
practical problem of the direct method is that it is typically an order of magnitude more 
demanding in computation than indirect method for a similar statistical accuracy. This 
is because for the direct method, many types (momenta, masses and sources) of quark 
propagators must be solved and a larger temporal lattice size is needed to accommodate the 
three point functions. 
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C. Lattice calculation setup 



We give here brief details of the aspects of the lattice calculation which are common to 
the following sections. More detail will be provided in each section. 



We use the domain-wall fermion (DWF) 



27 
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action for quarks. At the expense of 



an additional fifth dimension, DWF formulati on p reserves the flavor and chiral symmetries 
of continuum QCD at finite lattice spacing [3(|. There are two main reasons to use a 
formulation with good chiral properties: 

1. Because of chiral symmetry, mixing between operators in different chiral multiplets is 
prohibited. As will be explained subsequently, for the particular case of interest this 
implies that the operators are renormalized multiplicatively, as in continuum QCD. 
Hence, simple and clean handling of the operator renormalization is possible. 

2. In lattice gauge simulations, one of the most important sources of systematic error 
is due to finite lattice spacing, a. Chiral symmetry disallows 0(a) scaling violations 
for both on- and off-shell Greens functions, which will participate in our estimate 
of the nucleon decay matrix elements. As such, it suggest a mild dependence of 
any observables on a. Indeed, DWF has shown good scaling behavior in various 
hadronic quantities H y, h I34I An important consequence of the fact that 



the propagator is off-shell improved is that non-perturbative renormalization becomes 
much simpler allowing for the possibility of significant reduction in systematic error. 

For a fifth dimension of finite extent, there is still some explicit breaking of chiral symmetry 
by the DWF action. However, it may be hoped that such breaking is small enough for 
computationally reasonable extents of the fifth dimension that it may be either ignored, or 
treated as a small correction. Indeed, as it will be shown later, mixing of the operators is 
absent to a high degree of accuracy with a finite, affordable size of the fifth dimension. 

A convenient quantity to parameterize the size of the explicit chiral symmetry breaking 
is the residual mass, m res , which acts as an additive renormalization to the fermion mass 
m f , and is defined by considering the Ward-Takahashi (WT) identity 

Wx)J b 5 (0)) = 2m f (JZ{x)J h M) + 2< J 5 » J 5 b (0)> + ^V 5 b (0)>, (25) 
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(t » 1), (26) 



where A^ is the flavor-non-singlet axial current defined using all bulk 5d fermion field and 
is a point-split bilinear operator in 4d sense. J5 is the pseudoscalar field constructed from 
the 4d quarks which are defined from the 5d domain-wall fermion field by taking the values 
at both walls, s = 0, L s — 1 (L s is the size of the fifth dimension). J 5? is similar to J5, 
but, is defined with the fermion field located at the mid-points (s = L s /2 — l,L s /2) of the 
fifth dimension. The definition of m res is such that the WT identity takes the same form 
as that in the continuum, with a shifted mass of rrif + m res \t4\. Hence, A® is often called 
the "conserved axial current". Note that, by the WT identity, at mj = — m res the pion 
mass must vanish, m 2 — > as (rrif + m rcs ) — * 0. For further details of our conventions and 
notation see [s^ . 

It has been demonstrated in quenched calculations that the chiral properties of DWF 



are improved on configuratio ns g enerated by improved gauge actions 



34j,|36(. In particular 



the DBW2 gauge action |37l . l3q is superior in the smallness of m res at a given L s 
use this DBW2 gauge action for both quenched and dynamical fermion simulations. The 
drawback of using the DBW2 action is that the sampling of different topological sectors 
becomes harder for finer lattices, which has already been observed for a ~ 0.1 fm j^. The 
main quenched calculation in this study uses a coarser lattice, with a = 0.15 fm, where the 
sampling problem is absent. We also use an ensemble with a = 0.1 fm, where we are able 
to overcome the sampling problem by setting larger separation in Monte Carlo time. 

Because of the smallness of the explicit chiral symmetry breaking, A a „ can be treated 
as a (partially) conserved current to a good precision. As such, one can calculate the 
renormalization of the naive local (non-conserved) axial vector current, A®, by taking the 
ratio 



Here, due to the point split nature of A a Ax), a linear combination of the displaced operators 
should be used to get rid of O(a) error and to reduce 0(a 2 ) error [25J. Precise details of 

nn 

this technique can be found in |32L I34j . The axial current renormalization calculated in this 
way is used as a building block of our non-perturbative renormalization of the nucleon decay 
operators. 

Table0summarizes the simulation parameters. Our main effort is devoted to the quenched 
simulation (Nf = 0) at a = 0.15 fm, where we perform both direct and indirect measurements 
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TABLE I: Domain-wall fermion simulation parameters for quenched (Nf = 0) and unquenched 
(Nf = 2) runs with DBW2 gauge action. Lattice spacing a indicates the approximate value with 
the p mass input at the chiral limit. Unrenormalized, approximate value of m rGS is presented 
in physical unit. More detailed value of them will be shown in the later sections. "#configs" 
shows the number of configurations analyzed in either matrix element (ME) or non-perturbative 
renormalization (NPR) calculation. Unquenched simulation takes three set of degenerate masses. 



N f 


a [fm] 


6/g 2 Ll x L T L s M 5 


m res [MeV] #configs(ME) #configs(NPR) 





0.15 


0.87 16 3 x 32 1 2 1.8 


1.3 


100 


51 





0.1 


1.04 16 3 x 32 16 1.7 


0.04 


400 


55 


2 


0.12 


0.8 16 3 x 32 12 1.8 


2.3 


94x3 


37-47 



of the nucleon decay matrix elements. Finer lattice with Nf = will be used to discuss the 
finite lattice spacing and volume effects. An investigation using dynamical fermion for the 
indirect calculation of nucleon decay matrix element has also been performed to evaluate 
the quenching effect. 



III. OPERATOR RENORMALIZATION 



In order to relate matrix elements obtained with the bare lattice operator to the con- 
tinuum counterpart in a given renormalization scheme, one needs a well prescribed way to 
renormalize the lattice operator. In the literature perturbation theory has been used for the 
renormalization of the nucleon decay operators on the lattice. However, lattice perturbation 
theory suffers from bad conver genc e, primarily due to the tadpole contribution. Mean field 
improved perturbation theory [39] works much better. However, when compared to non- 
perturbative extractions there are some operators for which it is inaccurate. Moreover, the 
definition is ambiguous. To be precise: there are several ways to select the mean field factor, 
and the resulting renormalization factor sometimes depends on the choice, and it is indeed 
the case for the nucleon decay operator for DWF [3] . These problems are naively expected 
to disappear in the continuum limit. However, it is preferable to avoid such an ambiguity. 

Non-perturbative renormalization (NPR) solves these problems. In this work we employ 
the non-perturbative, MOM-scheme, renormalization technique of the Rome-Southampton 
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421 in 



group j^jj], which has previously been successfully used in conjunction with DWF 
the context of the renormalization of the flavor non-singlet quark bilinear operators and four- 
quark operators relevant to Kaon physics. Since the MOM- scheme can be applied to any 
regularization, it is also referred to as the regularization independent (RI/MOM) scheme. 
Our approach in this work is to use NPR on the lattice to extract the renormalization 
factors defined at some scale in the MOM-scheme, and then match the MOM-scheme to the 



MS'-scheme, which is more commonly used for the calculation of Wilson coefficients, using 
continuum perturbation theory. The renormalized operators are regularization independent 
up to discretization error for the lattice calculation, e.g. 0(a) for Wilson, 0(a 2 ) for DWF, 
and up to the truncation error in continuum perturbation theory. Since we are relying - 
in part - on continuum perturbation theory, we must work at a momentum scale for which 
it is applicable, and so we need this renormalization scale to be much greater than Aqcd- 
This may cause a problem because the lattice discretization error grows as the momentum 
becomes larger. This is the well-known window problem: is it possible to have a region of 
momentum such that Aqcd <C \p\ <C ft/a. We will give an estimate of the systematic error 
due to the window problem later. 

The anomalous dimension of the nucleon decay operator has been calculated up to two 
loops (43I using naive dimensional regularization (NDR) in QCD. We use this result both in 
the scheme-matching calculation, and to factorize the proper scale dependence of the NPR- 
MOM renormalization factor. In general, scheme dependence appears at the next-to-leading 
oder (NLO), and so the one-loop matching factor which relates the MOM scheme to MS, 
NDR is needed for the complete NLO treatment of the operator renormalization. The result 
for this will be presented in Appendix 



A. Operator mixing 

Given the good chiral symmetry of DWF, the mixing of the operators with different 
chirality is expected to be suppressed. However, it is instructive to first enumerate the 
allowed mixings if chiral symmetry is not assumed. Since this discussion uses only the 
rotational, parity, and the vector flavor symmetry of lattice QCD, it also gives the operator 
mixing structure for Wilson fermions. 

It is convenient to introduce the operator basis which mimics those commonly used for 
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TABLE II: Classification of the nucleoli decay three quark operator by parity (V) and switch- 
ing (5) (u <-> d). 





s- 


5+ 


v- 


SS PP 


AA 


VV TT 


v + 


SP PS 


-AV 


-VA TT 



the four- Fermi operators in the weak effective Hamiltonian. All the operators can be written 
in the form 

= e ljk (u iT CTd ] )T's k . (28) 

which should be Lorentz spinor, thus all suffixes other than single spin index must be con- 
tracted. Here, again u, d, and s are not necessarily labeling the real flavors. With a notation: 
S = 1, P = 75, V = 7 M , A = 7^75, T = = \{l^lv}, f = 7 5 a^, we have IT' = SS, 
PP, VV, AA, TT for the negative parity (V~) operators, and SP, PS, VA, AV, TT for 
the positive parity (V + ) operators [Zf|. There is another global symmetry which is useful 
in classifying these operators: switching (S) u and d is a symmetry of the Lagrangian if 
they are degenerate in mass. Under a switching transformation, an operator comes back to 
itself with possible change of sign depending on V that connects the spin indices of u and d. 
Recalling (CT) T = -CT for V = S, P, A (S~) and (CT) T = +CT for T = V, T (S+), we 
have four different operator groups as shown in Table |H] Operators in different blocks do 
not mix each other. 

These operators have the following properties: 

1. There is a trivial relation between operators with different parity in the same column 
in Table HP 0(V + ) = 750('P~). This means that there is a one to one mapping 
between the parity negative and parity positive operators such that the renormalization 
matrices are identical. 

2. The five operators in 0(V~) form a complete set of operators made of u, d and s with 
any ordering. This follows from the fact that any such operator can be rewritten, by 
Fierz transformation, as a linear combination of the operators 0^ s . 
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3. As our target operator is IT' = Pr/lPl (Eqs. (J3J), (jTJ)), we may neglect S + sector from 
possible mixing candidates. Then for each parity (chirality), we need to consider only 
three operators for mixing. 

4. Operators of the type udu are renormalized in the same way as uds. A simple way 
to see this is to note that the calculation of these renormalization factor using the 
Rome-Southampton NPR method is identical. This will be shown below. 



From now on, we can concentrate on the S sector. Our renormalization convention is, 

where a and b stand for possible IT" = SS, PP, AA for V~ . Zf^ D is a 3 x 3 renormalization 
matrix. As mentioned above, an identical matrix applies for V + operators, SP, PS, and 
AV. The chirality basis, which is more convenient to match the lattice operators with those 
used in the principal matrix elements, is made of following three operators, 

LL = ^(SS + PP) - ^(SP + PS), (30) 
RL = ^(SS - PP) - ^(SP - PS), (31) 
A(LV) = l -AA-\{-AV). (32) 

The renormalization matrix transforms under this basis change as, 

rvchiral rp yparity rp— 1 /qq\ 

ND — 1 ^ND 1 5 \ 66 ) 



t 1/4 1/4 » 



T 



(34) 



1/4 -1/4 
v 1/2 J 

where Z$g al is for the basis operators LL, RL, A(LV), while Z p ^ ty for SS, PP, AA. If 
there is no explicit chiral symmetry breaking by the action used, Z^ al is a diagonal matrix. 

Taking the above into consideration, we may contrast the situation when using Wilson 
and Domain Wall fermions: In the Wilson fermion case, lattice operator, for example LL 



operator, is renormalized at one loop 
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as 



Of/M = (l + ^[4MH + A]) &H + ^-(C RL O l $ + C A{LV) 0% V) ), (35) 



1(3 



where, A, Crl and Ca(lv) are scale independent constants. The discussion here has shown 
that these three operators in rhs are all that would appear even to all order of a s . In the 
DWF case, the mixing between operators with different chirality (off-diagonal elements of 
^ND al ) i s highly suppressed. This follows from the discussion which makes use of the low 



energy effective theory of DWF 4l|. Applying the same procedure, it can be simply shown 
that the off-diagonal elements are suppressed by a factor of m^ es (< O(10 -6 ) for all our 
parameters) . 



B. NPR formulation 



4l| we first calculate 



To calculate the renormalization factors using the MOM-scheme 
the Greens function of the particular operator in question with three external quark states 
in Landau gauge, 

G a (x ,x 1 ,x 2 ,x 3 ) = (Ol ds (x )u(x 1 )d(x 2 )s(x 3 )). (36) 

We set xq = 0. A Fourier transformation is then performed on the three external quark legs 
with the same momentum p, which are then amputated to obtain the vertex function, 

AV) = F.T. G°(0, 0:1,0:2,0:3)1^. (37) 

Writing the tensor indices explicitly, the renormalization condition of the MOM scheme 
reads, 

where Z q is the quark wave function renormalization; k are color indices, a, (3 and 7, 5 
are Dirac indices associated with T, T' respectively. 

The projection matrix P a is chosen such that Eq. (|38|) holds for the free field case with 
Z g = 1, Z a N b D = 5 ab . Then, 

TsC-y^s 7 (40) 

d^c-y a ( l5ltl ) 5 \ (4i) 

To simplify notation, we define the matrix M as, 

M ab = A^- fe a p 1& ■ P b j k p a Sj , (42) 
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pSS _ 






96 


ppp _ 






96 6 


pAA = 


1 




384 



which is equal to Z^ 2 (Z~^ l D ) up to systematic errors. 

The treatment of Z q needs care j^], as its definition naturally involves the derivative 
with respect to the momentum, and the lattice momentum cannot be continuous. Here we 
exploit the accurate determination of Za from the hadronic matrix element (Eq. (|2Tj) ). The 
bilinear vertex function of the local axial current calculated in the MOM-scheme yields 
Z q /ZA- Thus, the ratio A 3 / 2 /M aa will give Znd/Z^ 2 . Znd is calculated with these two 
measurements without directly dealing with Z q . 

As Za has no scale dependence in the continuum, it must not have scale dependence 
except for that brought by the discretization error ( which starts at 0(p 2 a 2 )) at the finite 
lattice spacing. The ratio Znd/Z^ 2 , then, has the same scale dependence as the nucleon 
decay operator up to 0(p 2 a 2 ) scaling violations. 



C. Scheme matching and RG running 

As it is clear in the above discussion, we need to know the scale dependence of the 
renormalized nucleon decay operator to separate it out from the potential lattice artifacts. 
Our goal is to quote values for the matrix elements of interest in the MS, NDR scheme at 
some scale, 

O^ifM) = u^s^^)o lat \ (43) 

where U MS ^ latt (fi) is the renormalization factor needed. We are using a two-step renormal- 
ization procedure: first renormalize with the MOM-scheme at scale p, and then match with 
the MS scheme at \i. This leads to the equation 

UMS^iau^ = U^; P )^^-Z^( P ). (44) 

The Z^ t and Zf£° are MOM-scheme factors calculated using continuum perturbation 
theory and NPR on the lattice respectively. Z MS (p) is the continuum MS renormalization 
factor, U MS (fi;p) is the renormalization group evolution factor from scale p to ji in MS. 
z JTs r ^ i 'z^° M \p) is the matching factor from MS to MOM at scale p. 

The anomalous dimension of the nucleon decay operator, which enters U MS (n;p), has 
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been calculated up to two loops in MS, NDR scheme 



7 
7o 



-4, 7l = - 



An 



The anomalous dimension reads, 

(45) 



14 4 ,r 

3 9 / 



4A 



where A = for LL, —10/3 for RL operator, and Nf is the number of active flavors 
The MS evolution factor reads 

« s (/i)l 7o/2fl) L . ( 7i A7o\ a a (ji) - a,(p) 



O 



«*(p) ' 

ii - In f , & = 102 - 



7i 
2 /3o 
22 
3~ 



2$ 



47T 



(46) 

(47) 
(48) 



The matching factor is calculated to one loop in continuum perturbation theory. The 
MOM-scheme calculation should be performed with the same kinematics and in the same 
gauge as that used on the lattice. Setting the momenta for the three external quark fields 
to be equal, and setting the mass to zero, the matching factor is obtained as, 



rMS 



gMOM 



4:71 



433 1123, 317, 

n 2 + £ In 2 

180 90 s V 180 90 



(49) 



where £ is the gauge parameter and £ = (Landau gauge) will be used. See appendixlUlfor 
the derivation. To match the NPR-MOM scheme with MS with this formula, we need to 
take the chiral limit of the massive simulation data. As we will see, this can be done very 
precisely, as this mass dependence is extremely mild in the NPR data. 



We use the two- 
for quenched QCD 



oop running coupling a s (/j) with A^g^ obtained by Alpha collaboration 



44j . A-^ro = 0.602(48), wher e rn is the Sommer parameter defined with 



MS 



the static quark potential V(r) as r 



2dV 
dr 



1.65 [45]. The approximate value is r = 0.5 



fm from the potential models. As we set the scale by using the p meson, we use our 
measurements oiro/a and m p a and combine with the experiment m p = 0.77 GeV to get the 
appropriate A^,. 



D. Results of the NPR 



We present here the results of the NPR of nucleon decay operators for quenched calcu- 
lation, on configurations generated with the DBW2 gauge action at a = 0.15 fm (see, Table 
P). The NPR study employs four quark masses mf — 0.025, 0.04, 0.055, 0.07, where the 
largest roughly corresponds to the strange quark mass. 

19 



Figure ^ shows the SS and PP projections of the SS operator, M ss,ss and M ss,pp 
as a function of lattice momentum squared for all quark masses. Note that mass depen- 
dence is negligible. Taking the chiral limit (mj- -* — m res ) using a linear extrapolation and 




ss,ss 



SS,PP 
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m f =0.04 
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m f =0.055 
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m f =0.07 
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FIG. 1: SS and PP projections of the SS operator, M ss,ss and M ' PP as a function of lattice 
momentum squared for each quark mass. 

rearranging to the chirality basis, one obtains Fig. |2 for all the elements of M. Most of 
the off-diagonal elements are consistent with zero for (pa) 2 > 1.2 within 2cr, while others 
still remain within 1% of the diagonal elements, and are thus negligible for our extraction 
|zs|. As a result, nucleon decay operator Q R / LL = e^ k (u lT CPR/L<V)PLS k is renormalized 
multiplicatively for our domain-wall fermion simulation. 

The next step is to obtain the total renormalization factor to relate the lattice operator 
to the MS, NDR operator, for which we need the value of Z q . As mentioned previously, we 
extract this value by calculating using the Rome-Southampton technique, and Za 

from hadronic correlators. For the former, we use the average of vertex function of the local 
axial vector and vector current operators. The renormalization constants for these operators 
are equal in a theory in which chiral symmetry is only softly broken. This equality should 
also hold for the vertex functions at high energies. At low energies they can differ due to the 
spontaneous breaking of chiral symmetry. Figure 01 shows the average and difference of the 
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diagonal 




O LL,LL 
□ RL,LL 
O A(LV),LL 
A LL,RL 

RL,RL 
V A(LV),RL 

LL,A(LV) 
+ RL,A(LV) 

A(LV),A(LV) 



FIG. 2: Mixing matrix M a,b in the chirality basis at the chiral limit rrif —> —m ies . 

vertex function for A and V. The non-zero difference at the value of the momentum which 
are accessible to our lattice calculations indicates that the window where the RI/MOM 
NPR can be safely applied is closing for the lattice spacing being used (a = 0.15 fm). This 
difference may be taken as a measure of the systematic error of the renormalization constant 
arising from the closing of the window. One sees up to 1.5% effect difference for p 2 > 1.2, 
which can be enhanced by the extrapolation, (pa) 2 — > to 2%. We may estimate the 
systematic error of the nucleon decay renormalization constants, from this source, as 3% 
considering the dimension of the operator. 

Using average of the axial and vector vertex function, combined with the vertex function 
of the nucleon decay operator, the diagonal elements of Zfj£ M / Z^[ 2 can be calculated. This 
is shown in Fig. |1] for O ll versus the renormalization scale squared. Using the results of 
the previously discussed matching and running calculations the MS renormalization factor 
(Eq. (|44jl ) at fixed scale \x = 1/a, which should, thus be scale invariant, is shown as squares. 
This is, again, plotted versus the square of the scale at which the lattice, MOM-scheme, 
renormalization calculation was performed. We identify the remaining momentum depen- 
dence as 0(p 2 a 2 ) discretization error. To extract the value at (pa) 2 = 0, linear function in 
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FIG. 3: Average and difference of the axial vector and vector vertex functions. 

(pa) 2 is used to fit it in the region 1.2 < (pa) 2 < 2.5, where the non-perturbative effect is 
expected to be small 3%) and higher order effect is negligible. Combining with the axial 
current renormalization ( Table liTT|) obtained with Eq. (|27|). and running from // = 1/a to 2 
GeV by Eq. (gTj), we get 



U 



MS<-latt 



(2GeV) 



0.751 (13) (45) for O ll , 
0.755(15) (45) for O rl , 



(50) 



where the value for O rl has similarly been calculated. The first brackets show statistical 
error. The systematic error in the second brackets involves two parts: one is from the 
window problem, which we already have estimated as 3%. The other is the perturbation 
theory error arising from truncating the higher order terms, which are of the order of a 2 . We 
take a 2 Av) — 0.05 at the smallest matching momentum (pa) 2 = 1.2 as the relative systematic 
error [7J|. Thus, the total systematic error of the NPR is 6%. 

We apply the same procedure to the a = 0.1 fm DBW2 lattices and get 



U 



MS^latt 



(2GeV) = 0.805(9)(32) 



(51) 
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FIG. 4: Renormalization factor of nucleon decay operator (Oil) normalized by the axial cur- 
rent renormalization constant Za, for a) (circle) MOM scheme Zf^ M (p)/Z^ 2 as a function of 
renormalization scale p 2 and for b) (square) MS scheme jj MS ^ latt ^ = l/a)/Z^ 2 (Eq. (|43jl ) as a 
function of MOM^ MS matching scale p 2 . 

for both O ll and O rl . The systematic error from the window problem is negligible. The 
perturbation theory error is 4%, which is counted as the total systematic error for the 
renormalization factor. 



IV. MATRIX ELEMENTS 

In this section the details of the matrix element calculation in the quenched approximation 
are given. 

A. Parameters and Lattice Scale 

The matrix elements and related hadron spectrum are calculated on the quenched DBW2 
configurations with lattice spacing a = 0.15 fm (see Table|l}. The 100 configurations used are 
separated by 200 iterations of four overrelaxation and one heatbath steps. The quark masses 
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are rrif = 0.02, 0.04, 0.06, 0.08. The strange quark mass point is approximately between the 
largest two masses used. We use quark propagators with anti-periodic boundary condition 
in temporal direction to measure Za, fn r es- For all the other quantities we average two 
quark propagators: one with periodic, the other with anti-periodic boundary condition in 
temporal direction. In this way we effectively double the temporal size to 64, for which we 
may safely neglect the effects of the hadrons moving around the temporal boundary in the 
three-point function measurements. 

Table UTTl shows the n, p, and nucleon masses, m res , local axial current renormalization Za, 
and pion decay constant f„ from the pseudoscalar two point function, all calculated with the 
degenerate quarks. Non-gaugefixed wall sources are used for f T . The other quantities 
use a quark source fixed into Coulomb gauge and uniformly distributed in a spatial cubic 
box, whose size is fixed to be 8 3 to optimize the nucleon signal. The hadron masses are 
obtained by the standard two-parameter correlated fit to the two point functions. m res and 
Za are calculated by the ratio method, Eqs. (}26|) . (|2*7j) . m res extrapolated to rrif = is used 
as the chiral limit point rrif — > — m rcs . In the table, chiral limit values from a linear fit in 
rrif are also shown. For pion mass, m\ is also fit to linear function of rrif . The lower bound 
of the fitting range is always the smallest mass, rrif = 0.02. For all measurements except 
for the p meson mass we exclude largest mass rrif = 0.08 from the chiral fit to stay in the 
region of good linearity. Taking the experimental p mass input (0.77 GeV) and comparing 
to the chiral limit the lattice cut off is determined, 

a -1 = 1.312(27) GeV. (52) 

If we took = 0.131 GeV as an input, we would obtain a^ 1 = 1.294(13) GeV which is 
consistent with the p input. For the nucleon mass = 0.938 GeV as an input, a~ l = 
1.249(32) GeV. This scale ambiguity (5%) between p and nucleon input, likely predominated 
by combination of extrapolation and quenching errors, will be used as the systematic error 
of the matrix elements in physical units. Using a = 0.150(3) fm from p input, the spatial 
lattice size is obtained as L a = 2.4 fm, which should be large enough for nucleon with 
our mass range j^, and can accommodate momentum small enough to directly reach the 
physical kinematics region of the nucleon to pion decay. 

The simple linear extrapolation of the pion mass to the chiral limit gives a non-zero 
value. Rather than evidence of numerically significant chiral symmetry breaking which is 
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TABLE III: The residual mass m res , axial current renormalization constant Za, hadron masses 
from the box source propagator, pion decay constant f n from the pseudoscalar two point function 
with gauge unfixed wall source propagator for 6/g 2 = 0.87 (a = 0.15 fm) quenched domain-wall 
fermions. The nucleon decay low energy parameters a and (3 with bare operator are shown, as 
well. Chiral extrapolation by uncorrelated linear fit in mj within the "range" of m/ are performed 
to get the "chiral limit", which is defined as rrif — > —m res , except for mj — > to get the chiral 
limit of m res itself. For m n , squared values are fitted and "chiral limit" shows the square root of 
extrapolated m 2 . All dimension-full quantities are in lattice units. Errors are by jackknife. 



m f 


m res (xl0- 3 ) 


Z A m p m N 


a(xl0~ 3 ) 


/3(xl0" 3 ) 


0.02 


1.171(29) 


0.78405(38) 0.3060(25) 0.637(16) 0.858(17) 0.1107(9) 


-6.85(66) 


7.56(76) 


0.04 


1.076(25) 


0.78770(28) 0.4172(22) 0.6866(69) 0.9698(83) 0.1205(9) 


-7.54(46) 


7.61(44) 


0.06 


1.005(22) 


0.79173(23) 0.5066(20) 0.7352(42) 1.0722(62) 0.1290(9) 


-8.28(42) 


8.27(39) 


0.08 


0.955(21) 


0.79609(21) 0.5849(18) 0.7829(32) 1.1651(52) 0.1368(9) 


-9.05(42) 


9.02(40) 


chiral limit 


1.248(34) 


0.77983(46) 0.083(11) 0.587(12) 0.751(19) 0.1012(10) 


-6.03(71) 


6.57(76) 


range 
X 2 /dof 


0.02-0.06 
0.14(44)/1 


0.02-0.06 0.02-0.06 0.02-0.08 0.02-0.06 0.02-0.06 0.02-0.08 0.02-0.08 
0.29(18)/1 0.22(15)/1 0.01(16)/2 0.13(29)/1 0.31(19)/1 0.005(51)/2 0.49(60)/2 



not already taken into account by the residual mass, we take this as a sign of non-analyticity 
at small mass. In fact, fitting using a formula including quenched chiral logarithms suggested 
by the quenched chiral perturbation theory (Q%PT) leads to a fit with reasonable \ 2 under 
constraint that m n vanishes at m/ = — m rcs . We use m/ = — m ICS as the chiral limit of 
DWF, thus as the physical u, d masses in our approximation, throughout the paper. Apart 
from the pion in the chiral limit, we use leading linear dependence on the quark mass to 
interpolate or extrapolate to the physical points for all the quantities. 

The strange mass point mi is obtained as m^mf ,mf) = 2m 2 K using the linear fit 
results to interpolate, where m K = 0.497 GeV, 

m f = 0.0675(30). (53) 

This will be used to interpolate the mass of s in the kaon for the form factors. 
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B. Hadronic matrix elements 



The low energy parameters a and (3 (Eqs. (|20|). (|21|l) of the nucleon decay chiral La- 
grangian are calculated through the ratio of the two point functions, 

A quark propagator with a cubic box source of 8 3 volume at t — to is used for both-two 
point functions in numerator and denominator. The factor \[Z~ V is the overlap of the proton 
interpolating field to the normalized proton state, 



(0\J p (0,0)\p) = ^Z p u p , (55) 

for the local proton interpolating field J p (x) and the proton spinor with the standard rela- 
tivistic normalization (see appendix El Eq. (jA2|0 . Z p could be measured via the amplitude 
of the single exponential fit to the point-point proton two point function. However, this 
fit is often problematic. Instead we estimate it through the ratio of the proton two-point 
functions with point and box source propagators, both with point sink. Asymptotic value of 
the ratio gives the ratio of the amplitudes with the point-pion and box-point propagators. 
Given the amplitude of the box-point propagator, which is more accurate than that with 
the point-point, Z p is finally obtained. 

FigureElshows the ratio R a (t) (Eq. (J53j) ) for m/ = 0.04 and 0.08. By fitting to a constant, 
ce(rrif) is obtained. Figure EH shows the fitted a at finite m/, which is extrapolated with linear 
function in m/ to the chiral limit jsoj ]. Our estimate of the low energy parameters with the 
operators renormalized at \i = 2 GeV is 

-a = 0.0100(12)(14)(6) GeV 3 , (56) 
P = 0.0108(13)(15)(7) GeV 3 , (57) 

where the first bracket shows the statistical error from the the bare matrix element, second 
is systematic error from the scale ambiguity, and the third comes from the total error of the 
renormalization factor. 

The direct method amounts to calculating a ratio of the three- and two-point func- 
tions. For the proton to ir case, the ratio is 

E^e^-^iUx^h) R / LL (x,t) J p (t )) 



E,^^^-)(j 7r (x 1 ,t 1 )4(f,t))-E.(^^)^(^o)) v 
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FIG. 5: Ratio R a {t) (Eq. (|54j)) for the low energy parameter a at mt = 0.04 and 0.08 on quenched 
DBW2 configurations, shown against the operator position (t), where the proton source is located 
at to = 6. 

with Z n being the overlap of the pion interpolating field to the normalized pion state, 

(tt|4(0)|0) = V^. (59) 

Again, the proton interpolating field at t — to is made of quark fields distributed uniformly 
in a 8 3 box. The associated quark propagators are solved with Vfif — m\. In the three-point 
function, as depicted in Fig. Q two quarks from the nucleon source are annihilated by the 
operator at t — t. The other is a spectator, which is annihilated by the pion interpolating field 
at t — t\. To interpolate the pion with momentum, a quark propagator is solved sequentially 
with m/ = m>2 and with the source equated with the spectator quark propagator at t = t\ 
with the momentum projection e ip ' x , which acts as an injection of the momentum p. The 
resulting sequential quark propagator is finally contracted as the third quark of the operator 
at t — t. The momentum — p is injected to the operator. 

The pion two point function in the denominator is made with the two quark propagators 
constructed using non-gaugefixed wall sources at t = t\ and with mass mi or mi. One source 
is always a zero momentum wall, while the other is made with distribution of e ip ' x , where 
the momentum p should be matched to the momentum injected to the sequential quark 
propagator in the numerator. Combining the two propagators of this type, the pseudoscalar 
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FIG. 6: | a | and \j3\ obtained from the ratio R a /p(t) as functions of m/ on quenched DBW2 
configurations. Linear extrapolation is performed to get the values in the chiral limit (m/ — > m res ). 
Values are in lattice units and unrenormalized. 

operator at the t\ timeslice becomes a local operator after averaging over the gauge config- 
urations, as the non-local terms vanish by Elitzur's theorem |48j. The non-gaugefixed 
source works well for the pseudoscalar operator and is superior to the local source in sig- 
nal/noise ratio. At time t, the quark and the anti-quark are annihilated by the local pion 
field Jl(x) with momentum injection. 

We only work with masses satisfying m\ < m 2 , which is enough to extrapolate/interpolate 
to the pion ((mi,rri2) — > {—m rcs , — m rcs )) or kaon (— m rcs , m^)) physical point. For extrap- 
olation to physical kinematics, pL a /2n = (1,0,0), and (1,1,0) are used, where L a = 16 is 
the spatial size. Zero momentum p = (0, 0, 0) data is taken only for mi = m 2 points to 
discuss the soft pion limit in Appendix iBl 

The matrix element is extracted by fitting the plateau of the ratio as a function of the 
position of the operator t. Some care is needed to extract the relevant form factor Wq Q|. 
The nucleon interpolating field also interpolates the parity partner of nucleon, which have to 
be eliminated by the parity projector. Then naively taking the trace of the ratio, we obtain 

Tr |p L [Wo - i4W q ] {^Y^j } = W - iq^Wg. (60) 
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FIG. 7: Diagram of the quark line contraction of the nucleon decay (N — > PS) three-point function. 



The irrelevant form factor is calculated as 



Tr\P L [W Q -i4W q 



1+74 



(61) 



where qj = —pj, the injected spatial momentum with negative sign. Combined with the 
calculated energy transfer zg 4 = E p — E n , Wq is finally disentangled. 
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FIG. 8: Ratio Eq. (f3H|) for the relevant form factor Wo of (tt°\ (ud)RUi\p) with bare operator and 
in the lattice unit at mj = 0.06. 

Figure |Hl plots the WVpart of the ratio Eq. (fo"%|) as a function of t for (7i°\(ud) R UL\p) at 
rrif = 0.06. Fitting plateau in the range 10 < t < 15, which appears to be a reasonably good 
fitting range for all the processes and parameters, Wo(rrif, q 2 ) is obtained and plotted in Fig. El 
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FIG. 9: Relevant form factor Wo of {nP\{ud)iiUL\p) as a function of — q 2 . Circles and squares 
are the measured data with pion having the smallest and the second smallest momentum on the 
periodic lattice respectively. Larger symbols show the degenerate mass points (mi = 1712), while 
the smaller indicate non-degenerate (m-i 7^ m.2). The line shows the chiral limit obtained from the 
fit (Eq. 1)62(1 ). The solid diamond indicates physical kinematics point q 2 — > 0. The open diamond 
is from the indirect method. 

as a function of — q 2 . The results are given in lattice units and for the bare operator. Solid 
diamond shows the extrapolation to the physical point (mi — > — m res , ^2 — * ~ m rcs, Q 2 ~^ 0), 
using jsi| 

W = c + ci(mi + m rcs ) + c 2 (m 2 + m res ) + c 3 q 2 . (62) 

This formula is obtained using leading order x?T, expanded in q 2 and quark mass m q = 
771/ + rn TCS . Quark mass is naturally interpreted as m 2 PS , which has similar size as q 2 at the 
simulated points. The open diamonds are from the indirect method, Eq. (f2*U|) taking f n = 
0.131 GeV from experiment. Inputting /„. from our lattice measurement gives a consistent 
result. Since the operators are renormalized multiplicatively (see Sec. IIH)) the results will be 
unchanged after renormalization up to an overall factor. The indirect calculation estimates 
Wo to be larger by about a factor of 2 than the direct method. Similar disagreement is 
seen for the LL operator case, which is shown in Fig. This difference, though large, 
is not surprising. At the physical kinematics point, where the pion momentum is large 
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FIG. 10: Same as Fig.H but for (n°\(ud) L u L \p)- 

\p\ = m^jl ~ rrtKi one might expect LO t° lose its effectiveness. On the other hand, 
the soft pion limit should be described by LO exactly. The discussion on this point is 
given in Appendix [5] 

The same analysis has been carried out for the other principal matrix elements. All the 
results for PV ( m i ; m 2, Q 2 ) are summarized in Table I VI II The fit results with Eq. (J52j) are 
listed in Table El Figure ^2 plots the final results of Wq at the physical kinematics in 
physical units for the renormalized operator. The indirect method estimates the form factor 
larger than the direct method for most of the cases. Note that, for several processes, the 
difference of the two estimates are significant. 

Table HVl shows the relevant form factors for all the possible matrix elements. This is the 
main result of this paper. We note that the total error is dominated by the statistical error 
of the matrix element. 

The values of individual matrix elements are different from those obtained for a -1 ~ 2.3 
GeV with Wilson fermion jxjj, while ratios of the matrix elements are similar. We expect 
our DWF results to be closer to the continuum limit, as discussed below. 



0.15 
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TABLE IV: The relevant form factors Wo [GeV 2 ] for all the possible matrix elements (Eqs. 
(|14[0 . with the operator renormalized at \x = 2 GeV with MS, NDR. The first error is total error, 
which is obtained by adding quadrature the statistical error of bare Wo (quoted second error), 
systematic error from the scale, and the total error of the renormalization factor. Note again that 
in the present calculation the flavor SU(3) breaking effect of r\ is not taken into account. 



Matrix element 


RL or LR operator 
Wq [GeV 2 ] total error stat. error 


LL or RR operator 
Wq [GeV 2 ] total error stat. error 


(it \(ud)u\ P ), 


(ir°\{du)d\n) 


-0.060 


0.018 


0.017 


0.086 


0.022 


0.019 


(<K+\(ud)d\ P ), 


— (tt + \(du)u\n) 


-0.085 


0.026 


0.024 


0.122 


0.030 


0.027 


(K°\(us)u\p), 


-(K+\(ds)d\n) 


0.082 


0.018 


0.015 


0.050 


0.012 


0.011 


(K+\(us)d\p), 


-{K°\(ds)u\n) 


-0.029 


0.008 


0.008 


0.028 


0.008 


0.007 


(K+\(ud)s\p), 


-{K°\(du)s\n) 


-0.090 


0.020 


0.017 


0.106 


0.021 


0.017 


(K+\(ds)u\p), 


-{K°\(us)d\n) 


-0.053 


0.012 


0.010 


-0.078 


0.015 


0.013 


(r]\(ud)u\p), 


— (r}\ (du)d\n) 


0.017 


0.008 


0.007 


0.078 


0.020 


0.017 



C. Systematic errors 

The major sources of systematic error for this calculation are: 

1. finite lattice spacing a, 

2. finite system volume, 

3. chiral extrapolation, 

4. quenching effects. 

Among these, Qand 121 are relatively easier to address, and, as we see below, there are reasons 
to believe they are small. 

To evaluate the systematic errors ^ and El we compare quantities already shown with 
those obtained on the finer lattice spacing and smaller volume. We have used a Q/g 2 = 0.87 
lattice, where lattice spacing is a = 0.15 fm, 16 3 lattice volume corresponds to (2.4 fm) 3 
box. The finer lattice parameter is a Q/g 2 = 1.04, a = 0.1 fm, (1.6 fm) 3 box. 
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FIG. 11: Relevant form factors of all the principal matrix elements for nucleon decay with the 
quenched approximation. Operators are renormalized by MS, NDR scheme at fi = 2 GeV. The 
errors are statistical only. 

First, we compare the hadron masses in Table ITTT1 with those obtained on the finer 
lattice. Figure El shows the hadron masses normalized by ro j3| against the quark mass 
r Z m (mf + m res ) renormalized at fi = 2 GeV in MS, NDR. The quark mass renormalization 
factors are taken from ref. j^J, where only L s = 16 case is studied. The difference of Z m 
between L s = 12 which has been used in this study and L s = 16 should be negligible. The 
error shown in the figure is statistical only. 

At larger masses hadrons are compact, hence the volume effect should be small. The 
consistency of the hadron masses at larger quark mass suggests the finite lattice spacing 
error is negligible within our accuracy. Furthermore, the consistency is seen all the way 
down to the smallest mass of a = 0.15 fm. This suggests the volume effect is also negligible 
in the mass range studied. 

More direct evidence that the systematic errors Q and 121 are small is seen in Fig.^3 where 
the difference of the two low energy parameters \a — f3\ is shown in the same way as Fig. 
We have used the relation of the proton interpolating field (Eq. ()A9j) ) to a nucleon decay 
operator, 

J P = (o*L - o L ±) + (©a - (63) 
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FIG. 12: Hadron masses as a function of quark mass renormalized at fj, = 2 GeV with MS, NDR. 
Both axes are normalized by r$. Both lattices have 16 3 x 32 volume. Larger lattice spacing a = 0.15 
fm has been used for the nucleon decay matrix element. Volume and lattice spacing effects are 
negligible. 



hence, 



\J p \p) = {a- (3)u p . 



(64) 



Thus, a — (3 can be extracted from the ordinary nucleon two point function (J p (x)J p (0)) as 
well. Since the renormalization factors of the nucleon decay operators of different chirality 
are same within the error, the average of them may be used to renormalize J p to get a — 3. 
For the finer lattice results we have reused the a = 0.1 fm data for the spectrum study |3jJ. 
The same arguments as the hadron masses lead to the observation that the systematic errors 
from lattice spacing and the volume are negligible also for a — (3. From this, we can expect 
that those errors on individual a, j3 and further the relevant form factor Wq are negligible 
compared to the other errors. 

Systematic error from the chiral extrapolation may be caused by the linear extrapolation 
in quark mass, while in general, contribution of higher order terms, and in the special case of 
quenching (by QxPT), lower order terms, may not be negligible. Our linear extrapolations 
for matrix elements and hadron masses (for extracting the lattice scale) have been performed 
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FIG. 13: Difference of the nucleon-decay low-energy parameter — (a — (3) (taking the sign of a and 
(3 from the previously descibed extraction) as a function of quark mass. Normalization is same 
as Fig. 1121 Consistent results are obtained for the different lattice spacings and different spatial 
volumes. 

for the mass region where the resulting % 2 is small. This does not necessarily mean the 
leading linear dependence is sufficient. This can be the case when a term with lower power 
occurs at higher loops in the chiral expansion. An example of this is the p j^] and nucleon 
masses, for which m l J 2 terms appear at one loop. The correct quenched results may only 
be found when using a proper Qx?T formula with data sufficiently good in quality and in 
quantity to determine the fitting parameters. Only after that we can deal with the quench 
error by examining the difference to the experiment or by comparing to the unquenched 
results. Lack of the QxPT knowledge of the matrix elements as well as the limited quality 
and quantity for our data make it difficult to follow this scheme. 

We observed a scale inconsistency between the p and nucleon mass inputs when using 
the linear chiral extrapolation. This difference has been taken as the systematic error to the 
matrix elements due to the scale ambiguity. Since the finite volume and finite lattice spacing 
errors are negligible, we expect this scale error is dominated by the chiral extrapolation and 
quenching error. Whether this estimate is plausible or not will be examined by comparing 
with the results with dynamical fermion simulation, which will be discussed in the next 
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section. Note, however, the comparison is done only at a finite lattice spacing for the 
unquenched calculation. While we would also expect the discretization error to be small for 
the unquenched calculation, this should be checked in future studies. 



V. DYNAMICAL QUARK EFFECTS 

Quenched calculations have, in general, an uncontrolled systematic uncertainty. In the 
previous section, we estimated the systematic error due to quenching as approximately the 
size of scale ambiguity. This is clearly an unsatisfactory technique, and full QCD calculations 
are necessary. Ultimately this can be completed by unquenching the u, d and s quarks in the 
direct calculation of the form factors. As a first step toward this direction, we shall examine 
the low energy parameters a and (3 with dynamical u and d quarks, while still treating the 
s-quark in the quenched approximation. 



A. Description of the simulation 

We use the two-flavor dynamical DWF configurations described in |53j. These lattices 
were generated using the DBW2 gauge action with Q/g 2 = 0.8, on lattices of size 16 3 x 32, 
a fifth dimension of L s = 12, and a domain wall height of M5 = 1.8 (see Table HJ). Periodic 
boundary condition was imposed for all except the temporal direction of dynamical fermions 
and valence fermions for the spectrum and matrix elements, where anti-periodic boundary 
conditions were used. The dynamical quark masses are m dyn = 0.02, 0.03, 0.04, which 
approximately covers the range from the half strange to strange quark mass. All the analysis 
was carried out for the valence masses equal to the dynamical masses m^ a/ = m dyn . 

For the mass and matrix element calculation Coulomb gauge fixed wall sources were used. 
Measurements were performed twice per lattice: with the two source time slices separated 
by the half size of the temporal direction L T /2, which effectively doubled the statistics in 
comparison to the analysis done in Ref. j&J. 

The results of the mass measurements are summarized in Table El A linear chiral ex- 
trapolation was performed for all the quantities in the table using all three quark masses. 
Figure El shows p and nucleon mass and their chiral extrapolations. 
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TABLE V: m rcs , Za and hadron masses from the wall source propagator for two- flavor dynamical 
domain- wall fermions with G/g 2 = 0.8 DBW2 gauge action. Two quark propagators with different 
source position are analyzed for each configuration. The nucleon decay low energy parameters a 
and j3 for the bare operator in lattice units are also shown. We have performed a linear chiral 
extrapolation for all quantities using all dynamical masses m/ = 0.02, 0.03, 0.04. The chiral limit 
is defined as rrif — > — m res , except mt — > for m ros . For m^, squared values are fit and chiral limit 
shows the extrapolated m 2 . 

mf rre rcs (xl0~ 3 ) Za i% m p tun a(xl0~ 3 ) /3(xl0~ 3 ) 

0.02 1.360(26) 0.76035(27) 0.2916(15) 0.5474(50) 0.7631(88) -4.77(25) 4.74(25) 

0.03 1.357(21) 0.76187(23) 0.3563(13) 0.5991(51) 0.8383(87) -5.27(25) 5.13(25) 

0.04 1.336(24) 0.76323(21) 0.4084(22) 0.6325(52) 0.8980(92) -5.97(31) 5.95(32) 

chiral limit 1.388(55) 0.7573(6) -0.0029(25) 0.459(11) 0.621(20) -3.47(62) 3.43(61) 

X 2 /dof 0.10/1 0.08/1 0.55/1 2.11/1 0.51/1 0.11/1 0.42/1 




0.01 0.02 0.03 0.04 

dyn 



FIG. 14: p and nucleon masses as functions of m?™ = m 
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The lattice spacing determined in the chiral limit by inputting m p = 0.77 GeV is 

a- 1 = 1.678(40) GeV. (65) 

which is consistent with, but, has improved statistical error to, the previous estimate gained 
using a single quark propagator on each lattice a" 1 = 1.691(53) GeV For comparison, 
a^ 1 = 1.688(21) (^94) is obtained [53[ with r = 0.5 fm input, which is consistent with p 



input. The lattice spacing from the p input is a p = 0.1176(28) fm. Thus, the spatial size 
of these lattices is 1.88 fm. The spatial size divided by the pion Compton wave length at 



our lightest mass is LM n ~ 4.7. Systematic numerical studies 



hi 



55|, as well as theoretical 



calculations |56j, on the finite volume effects in the nucleon mass with two flavor Wilson 
fermions indicate just a few percent finite volume mass shift with similar parameters to our 
lightest point. By this we expect the volume effect on the low energy parameters are also 
small and negligible compared with the statistical error. 

The ratio m^jm p = 1.353(47) at the chiral limit is larger by 3a from physical value 1.218. 
While this discrepancy is of a size that could be easily ascribed to scaling violations, it should 
be noted that a similar size of discrepancy for the ratio has been observed in the continuum 
limit of the two-flavor Wilson fermion [5j| using the polynomial chiral extrapolation for 
the masses. This suggests that our discrepancy might persist towards the continuum limit. 



rer order chiral expansion for the 



There is a way to cure this problem by employing the hig 
nucleon mass with terms non-analytic in the quark mass 

B. Matrix elements 

For the non-perturbative renormalization of the operators we follow the same procedure 
as for the quenched calculation. The dynamical mass points mj = 0.02, 0.03, 0.04 are used 
to analyze for the MOM scheme renormalization. Linear chiral extrapolation in dynamical 
mass mj is carried out. The resulting renormalization factors (Eq. f!43|) ) which renormalize 
the lattice operator to give those in MS, NDR at \i — 2 GeV are 

. 0.731(28)(39) for C i;L , , 
U MS ^ latt (2GeV) = { (66) 

0.722(28) (39) for O r ' l , 



where the calculated by Alpha collaboration for the two flavor QCD 
0.62(4) (4) has been used. 
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The low energy parameters are also obtained in the same way as in the quenched case, 
and are shown in Table Figure shows a at the renormalization scale p = 2 GeV 
obtained on the dynamical configurations as a function of the pion mass squared with scale 
set by p in the chiral limit. The quenched results are shown for comparison. The dynamical 
result has stronger m n dependence than quenched. After rather long extrapolation to the 
chiral limit with linear function of quark mass we obtain the a and j3 parameters as 



a = 0.0118(21) GeV 3 , 
P = 0.0118(21) GeV 3 , 



(67) 
(68) 



where errors are only statistical. Comparing to these large errors, the errors of the renormal- 
ization constants are negligible. These values are consistent with those obtained in quenched 
approximation Eqs. (JHfij) - (J57j) . 

T 




rri 2 [GeV 2 ] 



FIG. 15: Comparison of the quenched and dynamical a as a function of pion mass squared. Scale 
is set by p input at the chiral limit. Operators are renormalized by MS, NDR scheme at p = 2 
GeV. 

In table I VII our results on the low energy parameters both for quenched and dynamical 
simulations are compared with those obtained in the literature, which include various QCD 
model calculations and the lattice QCD efforts by several groups. It should be remarked 
that by the efforts of recent calculations including this work substantially have decreased the 
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TABLE VI: Comparison of the low energy parameter of the nucleon decay chiral Lagrangian a 
and (3 among various QCD model calculation, lattice results in the literatures and the results from 
this work. In lattice QCD calculations, WF and DWF mean Wilson and domain-wall fermions. 
Our quenched results are shown with the total error consisting of statistical and systematic errors 
on the bare matrix elemnt, renormalization constant and scale. The unquenched errors are only 
statistical. 

H [GeV 3 ] \f3\ [GeV 3 ] 



Donoghue and Goldwich [10] 0.003 




Bag model 


Thomas and McKellar [L3J 


0.02 




Bag model 


Meljanac et al. [IT] 


0.004 




Bag model 


QCD model Ioffe [8] 


0.009 




Sum rule 


calculation Krasnikov et al. [12] 


0.003 




Sum rule 


Ioffe and Smilga [14] 


0.006 




Sum rule 


Tomozawa [9] 


0.006 




Quark model 


Brodsky et al. [15] 


0.03 






Hara et al. [16] 


0.03 




WF, a = 0.11 fm 


Bowler et al. [171 


0.013 


0.010 


WF, a = 0.22 fm 


Lattice QCD Gavela et al. [18] 


0.0056(8) 


~ a 


WF, a = 0.09 fm 


N f = JLQCD [19] 


0.015(1) 


0.014(1) 


WF, a = 0.09 fm 



CP-PACS & JLQCD [20] 0.0090(09) (tig) 0.0096(09) (t 6 20 ) WF, continuum limit 
This work 0.0100(19) 0.0108(21) DWF, a = 0.15 fm 

Lattice QCD 

This work 0.0118(21) 0.0118(21) DWF, a = 0.12 fm 

Nf = 2 

ambiguity of the low energy parameters a and (3. Their absolute values fall in the range of 
the approximate value 0.01 GeV 3 , which lies in the middle of the wide range of the various 
QCD model calculation result 0.003 - 0.03 GeV 3 . 
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VI. CONCLUSIONS 



Using domain-wall fermions and a non-perturbative renormalization we have calculated 
the low energy matrix elements of the nucleon decay whose operator would appear in the 
lowest order, low energy effective Hamiltonian of any fundamental high energy theory that 
breaks baryon number and respects the symmetries of the standard model at low energies. 

The low energy parameters a and j3 are important fundamental quantities of the SU(3) / 
baryon chiral Lagrangian with the baryon number violating interaction. In the quenched 
approximation they were estimated at our lattice spacing (a ~ 1.5 fm) as a = —0.0100(19), 
(3 = 0.0108(21) GeV 3 for the operator renormalized at fi = 2 GeV. We have examined the 
finite volume and finite lattice spacing effects on a — j3, which appeared to be negligible. 
These results are consistent with those using Wilson fermions in the continuum limit j^J. 
The dynamical quark effects on a and (3 were examined by unquenching u and d quarks. 
Their values at the chiral limit are consistent with those in the quenched approximation. 
All the results so far indicate \a\ ~ \(3\ ~ 0.01 GeV 3 , which lies in the middle of the various 
model calculations (0.003 - 0.03 GeV 3 ). 

However, these parameters are not quite useful for the physical decay process. The results 
of the nucleon decay form factors calculated with a and (5 showed clear deviation from the 
direct calculation for several processes. This is presumably due to the large energy of the 
pseudoscalar meson, and because of neglecting some leading quark mass dependent terms 
in the chiral Lagrangian. Instead, it is advised to use the results with the direct calculation 
fTable HV)) for the partial width of proton, Eq. (|18|) . The value of the form factor with 
the direct method tends to be smaller compared to the indirect method; thus the direct 
method, which is more accurate, tends to prolong the proton life-time. Indeed, for the case 
of p — > e + + 7r°, experimentally which is important especially for model independent analysis 
such like Ref. |60j], the difference in our central values (with appreciable error bars) for the 
form factors tend to approach factor of about two. It is clearly important to improve the 
precision of the direct method in the future calculations. 

The dynamical quark effects on the low energy parameters the mass dependence became 
stronger on unquenching. However, this comparison was performed at a single lattice spac- 
ing, and with only two - relatively heavy - dynamical quarks. A clear target for the future 
is to repeat this calculation with three flavors of dynamical quarks with masses near the 
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physical values and multiple lattice spacings. 
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APPENDIX A: FIELDS AND CONVENTION 



In this appendix, various conventions in our calculation are summarized. 
Let p and k be spatial momenta and % and j symbolically denoting the other discrete 
quantum numbers, the state normalization is given by 

{p,i\k,j) = (2 7 r) 3 2E{p,i)6 3 {k-p)6 i>j , (Al) 

where E(p, i) is the energy of the (p, i) state. 

Spin 1/2 wave function uat(/c, s) with momentum k, spin s, and mass m, that obeys the 
Dirac equation (—itfiiN = mux), has a relativistic normalization, 

TiN(k, s)uN(k, s') = 2m5 SiS '. (A2) 
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Hadron interpolation fields, Juad are given as 



J n 


V ^ 


(A3) 


J-K+ 


= d^u, 


(A4) 


J-ft- 


= u<y 5 d, 


(A5) 


J K 


= sj 5 d, 


(A6) 


J K + 


= s^u, 


(A7) 


J II 


= -j=(u^ 5 u + d^ 5 d - 2s7 5 s), 
V6 


(A8) 


V 


= e ijk (u iT C-f 5 d j )u k , 


(A9) 


J II 


= e ijk (u iT C^d j )d k . 


(A10) 



APPENDIX B: CHIRAL PERTURBATION 

We use the tree-level chiral perturbation theory for the nucleon decay matrix elements 
[lil l^ . The strong interaction chiral Lagrangian of SU(3)f octet baryons has coupling 
constants: D and F describing the coupling of the nucleon to the axial current, aj, fej, where 
% = 1, 2, dealing with the leading linear quark mass dependence. The nucleon decay specific 
term is proportional to a or (3 that we are calculating on the lattice. 

The low energy parameters D and F can be related to the axial charge associated with 
the baryon semi-leptonic beta decay. D + F = g^ ^ = 1.27 is the nucleon axial charge, 
while D-F = gf~ n) = 0.33 - 0.34 can be measured by S — > n + e + u e dj are 

determined by the mass difference of the baryons, and the size is 0(1). For nucleon decay 
matrix element, it enters through the baryon masses. While hi have direct influence to the 
nucleon decay amplitude, their values are not well-known. Although they can naturally be 
0(1), we set them to be zero in the following. 

Here is the summary of the value of parameters we employed, 



D = 0.8, (Bl) 

F = 0.47, (B2) 

/ = 0.131GeV, (B3) 

m N = 0.94GeV, (B4) 

m B = 1.15GeV, (B5) 
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where mg is the average baryon mass mg ~ Ms — M\. 

Using these constants and the approximations: m u ^ <C m s <C mjv,B, — q 2 <C m Ar,_B) the 
relevant form factors (PS\O\N) = Wq of all the principal matrix elements are obtained as 



7T 
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u, d) L u L 

u, s) R u L 
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\p)o 
\p)o 
\p)o 
\p)o 
\p)o 
\p)o 
\p)o 
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\p)o 
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V2f 

g 

a 
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m N 



en 



a 2D m N 
/ 3 m B ' 
(3 2Dm N 
f 3 ms' 

/( 1 + U 

(3 A / -D 



D — F) 



m B 
m N 
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f 

a 

1 



3 



F 
F 
F 



F 



;i + F-3F), 



m N 
m B 
m N 
m B 
m N 
m B 
m N 
m B 



VQf 



(3-F + 3F). 



(B6) 

(B7) 

(B8) 

(B9) 

(BIO) 

(Bll) 

(B12) 

(B13) 

(B14) 

(B15) 
(B16) 

(B17) 



Nucleon to pseudoscalar and nucleon to vacuum matrix elements are related in the zero 
momentum limit of the pseudoscalar by soft pion theorem, 



\im(7i k ;p,\0\N) = -- 



(0\[Q k 5 ,O]\N). 



(B18) 



Q k is the axial charge having the same £77(3) flavor content as the n h pion. It is one of 
all the pseudoscalar states: 7r 0,+ '~, K°' + , r], as we can consider an ideal situation, mass-less 
limit of all of them. 

By Eq. the matrix element in the soft pion limit is written in terms of the form 

factors as 

\im(PS; Pl ,\0 L \N} = P L [W + m N W q ]u N . (B19) 
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Following relations for (PS\ 0\N) sp = Wo + ra^W q are obtained by the soft pion theorem 
Eq. (|B18|) and also by p^ — > limit of the tree-level results of the chiral perturbation theory. 



(n°\(u,d) R u L \p) sp = (B20) 

(7r°\(u,d) L u L \p) sp = (B21) 

(K°\(u,s) R u L \p) sp = --, (B22) 

(K°\( U ,s) L u L \p} sp = (B23) 

(K + \(u,s) R d L \p) sp = 0, (B24) 

(K + \(u,s) L d L \p) sp = 0, (B25) 

(Y 

(K + \(u,d) R s L \p) sp = -, (B26) 

(K+\(u,d) L s L \p) sp = ^ (B27) 

(^ + |K S ) R n L |p) sp = -, (B28) 

(K+\(d,s) L u L \p) sp = (B29) 

a 

(r)\(u,d) R u L \p) sp = — — , (B30) 

3(3 

{r]\(u,d) L u L \p) sp = — =-. (B31) 



These exact relations no longer involve any ordinary low energy baryonic constants. 
Checking these relations with the lattice simulation provides a consistency test of the whole 
procedure. We have measured Wo—iq^Wq with zero pseudoscalar momentum and degenerate 
masses for Ihs. The chiral limit should be taken with 

W - iq^Wg = c + cim g + c h ^/m^. (B32) 

The square root of the quark mass (oc m,) has entered from q 2 dependence. Together with a 
and (3 (Eqs. (|56j). (|57j)) the difference A = Ihs — rhs is calculated and shown in Fig. 1161 Most 
of the processes are consistent with A = (note that these are highly correlated values). 
Even in the worst case the deviation is less than 2a. 
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soft pion limit, direct-indirect 



■<7t \(ud)„u,\p> 

o 

<7l \(ud) L u L \p> 
<K°\(us) R u L \p> 
<K°\(us) I u l \p> 
<K + \(us) R d ] \p> 
<K + \(us) L d L \p> 
<K + \(ud) R s L \p> 
<K + \(ud) L s L \p> 
<K + \(ds) R u L \p> 
<K + \(ds) L u L \p> 
<r\\(ud)„u \p> 
<r\\(ud) L u \p> 




-0.1 



0.05 

AGeV 2 



FIG. 16: The difference A of the form factors for direct and indirect method in the soft pion limit. 
It shows the consistency with the soft pion theorem. 




FIG. 17: The tree and three one-loop diagrams of the vertex function for e* jfc (TiP l rd : ')T' s k . 
APPENDIX C: PERTURB ATI VE MATCHING FACTOR 

Here we give a summary of the nucleon decay operator matching calculation. Through- 
out we make use of the Minkowski space path integral. We need to match to the lattice 
MOM scheme with the continuum MS scheme. Since the MOM scheme is regularization 
independent, we can use perturbation theory in both schemes. In the lattice MOM scheme, 
the same momentum p E is injected to all the three quark external lines, where (p E ) 2 > 0. 
The corresponding Minkowski momentum is (p M ) 2 = — (p E ) 2 < 0. Both the MS and MOM 
schemes are denned in the mass-less limit. We set all masses to zero for MOM scheme in 
perturbation theory. 
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The tree-level vertex function is 



A^ = e ^r®r', 



(Cl) 



where T and T' can be either 1 or 75 or their linear combinations. The corresponding 
Feynman diagram is shown in Fig. El (leftmost). The upper left fermion line is for the 
charge conjugated field, whose vertex for one gluon emission is —ig{t a ) T whereas the normal 
vertex is igt a . 

The sum of the three one-loop diagrams in Fig. El with the naive dimensional regular- 
ization (NDR) scheme reads 



A 



(i) 

ND 



(C2) 



where /1 is the arbitrary mass scale introduced by the dimensional regularization, e = (4 - 
D)/2, D is the regularized dimension, £ is the gauge parameter (£ = is the Landau gauge 
X F and \ L ~ F read 

(D + 2)(p 2 - k 2 )r®r 



x l-f 



[(p - k) 2 + ie] [{p + k) 2 + ie] [k 2 + ie] 

3r<g>r' 

[(p - A;) 2 + ie][(p + A;) 2 + ie] 

p 2 T ® r 

' [(p - A;) 2 + ie] [(p + A;) 2 + ie] [k 2 + ie] 

2(p • fc)r <g> r'jrf# 

' [(p - A;) 2 + ie] [(p + A:) 2 + ie] [P + ie} 2 ' 
Carrying out the momentum integral and identifying fi 2 with — p 2 , one finally finds 



(C3) 



(C4) 



A 



(i) 



a. 



^(f-H- (i -°®-W} Affi " (C5) 



where e 1 = 2/(4 — D) — 7^ + In 47r. The renormalization condition of the MOM scheme up 
to one loop reads 

A(0) _ ( yMOM\-3/2 7 MOM( a (0) , * (1) x / pR \ 

The MOM scheme quark wave function renormalization is given as (see for example 



7 MOM 



Air 



Then the nucleon decay operator renormalization factor is calculated as 



7 mom 

J ND 



47T 



, 433 1123, /587 317 

2e~ + — - In 2 + £ I — - — In 2 



180 



90 



180 90 



(C7) 



(C8) 
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The MS renormalization factor is given by taking only the term proportional to e 1 in 
the square bracket. The matching renormalization factor z M /Z has been shown in 
Eq. (Pf. 

APPENDIX D: NUMERICAL TABLES 

All the results for the relevant form factor Wq of the principal matrix elements are shown 
in Tables IVIIHIXI The fit result of the relevant form factors with Eq. (jo^j) is presented in 
Table H 
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Note that one can still save the theory by appropriate changes in the flavor structure 
Abbott and Wise reduced Weinberg's six operators to these four, using (anti) symmetric 
property under flavor exchange. 

We call all the qqql operators as dimension six. The naming is just from the naive counting 
of the field dimension and nothing to do with the origin of the operator. 
It requires to calculate disconnected diagram, which is challenging and demanding in the 
numerical simulation. 

From the numerical calculation where we set nucleon at rest, we observe Wq ~ —iqoW q . Thus, 
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with the existence of lepton spinor v c at left, v°(— i$W q ) = v c miW q ~ v c mi/qo ■ Wo. Since the 
qo is the order of the nucleon mass, —i(jjW q is smaller than Wq by an factor ~ mijm^. 

[73] As stated, direct calculation of the form factor does not take into account the flavor SU(3) 
breaking effect of r/. Thus, this 200% correction can be regarded as an estimate of the system- 
atic error of direct calculation of (r]\(ud)jiUL\p) ■ The correction for (j]\{u(1)lul\p) is 30%. 

[74] For the ratio of correlators given on the right-hand side of Eq. ()26|) to be represented by a 
single number, we are assuming that the effect of J§ q on the 4-dimensional physics is local, 
and we are neglecting effects of 0(a 2 ). 

[75] As with the extraction of the residual mass, for this equation to be true we are assuming 
locality and neglecting 0(a?) effects, and so there is some minimum value of t before which 
Za may not be extracted. 

[76] In this notation, the PT' operator consists of scalar diquark. 

[77] Typo in 2] corrected by private communication with Nihei. 

[78] It should be noted that we do not ascribe the non-vanishing values to be a consequence of 
explicit chiral symmetry breaking due to a finite fifth dimension, but a combination of low 
energy spontaneous chiral symmetry breaking effects and discretization errors. 

[79] Another way to estimate the truncation error is to guess from the difference between LO and 
NLO. It is about 15%, which can be seen in Fig. ^ The difference between NLO and NNLO 
is, then, estimated by taking a square as 2%. This is same order, but, smaller than 5% that 
we have employed. 

[80] One can take the ratio a/f^ ((3/ f n ) at finite mass with using measured and extrapolate 
it to the chiral limit instead of extrapolating a {(3) itself. This method has an advantage of 
making the mass dependence milder. With our measurement, however, it makes only up to 
2% difference in the extrapolated values and up to 10% difference in the statistical error. 
Apart from the way it is extrapolated, this ratio in the chiral limit might be more useful 
than the original, since a or (3 enters the form factors (Eqs. (|20j) . (|21j) . HB6j) - (|B17jO always 
through the ratio with /„-. Also, since the ratio has one dimension less than the original, the 
scaling violation is expected to be smaller. Thus, taking the continuum limit with this ratio 
has an advantage in the discretization where the scaling violation is larger than the present 
calculation. 

[81] JLQCD 19] included a q 4 term. With our data the effect of including such a term is to shift Co 
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within one third of the statistical error. While x 2 /dof of the uncorrelated fit increases slightly 
by eliminating q 4 term, the relative error of x 2 estimated by jackknife stays in a similar value. 



52 



TABLE VII: Relevant form factor Wo of the principal matrix elements in lattice unit with the 
unrenormalized operator for each mass and momentum calculated on the quenched a = 0.15 fm 
configurations. p\ = (1,0,0) • 2TxjL a is for the smallest, f>2 = (1,1,0) • 2ir/ 'L a is for the second 
smallest momentum. 



777-1 777,2 



0.02 0.02 
0.02 0.04 
0.02 0.06 
0.02 0.08 
0.04 0.04 
0.04 0.06 
0.04 0.08 
0.06 0.06 
0.06 0.08 
0.08 0.08 



0.02 0.02 
0.02 0.04 
0.02 0.06 
0.02 0.08 
0.04 0.04 
0.04 0.06 
0.04 0.08 
0.06 0.06 
0.06 0.08 
0.08 0.08 



RL operator LL operator 

P = Pl P = P2 P = Pl P=P2 



(<K°\(ud)u\p) 
-0.053(15) 0.051(44) 0.061(19) 0.057(72) 
-0.060(12) 0.009(21) 0.065(14) 0.055(35) 
-0.064(11) -0.011(15) 0.067(13) 0.059(24) 
-0.066(11) -0.023(13) 0.069(12) 0.060(18) 
-0.071(8) -0.036(10)0.075(8) 0.061(14) 
-0.073(7) -0.042(8) 0.075(7) 0.060(10) 
-0.073(6) -0.047(7) 0.074(7) 0.061(8) 
-0.082(6) -0.056(7) 0.082(6) 0.064(7) 
-0.081(5) -0.058(6) 0.080(5) 0.064(6) 
-0.088(5) -0.065(5) 0.086(5) 0.068(5) 



(K°\(us)u\p) 
0.054(13) 0.050(43) 0.035(14) 0.019(49) 
0.056(11) 0.035(22) 0.036(10) 0.016(24) 



0.057(11) 0.035(17)0.039(9 

0.059(10) 0.038(14)0.042(8 

0.065(7) 0.045(11)0.031(5 

0.066(7) 0.047(9) 0.034(4 

0.066(6) 0.050(8) 0.037(4 

0.069(5) 0.053(7) 0.033(3 

0.069(5) 0.055(6) 0.036(3 

0.072(4) 0.058(5) 0.036(2 



0.019(17) 

0.023(14) 

0.017(10) 

0.024(8) 

0.029(6) 

0.029(5) 

0.033(4) 

0.034(4) 
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TABLE VIII: Continued from Table Eg} 



RL operator LL operator 

P=Pl P = P2 P = Pl P = P2 



(K+\(us)d\p) 



-0.020(9) 


0.013(29) 0.017(10) 0.017(35) 


-0.020(7) 


0.002(13)0.018(7) 


0.015(16) 


-0.020(6) 


-0.001(9) 


0.019(6) 


0.016(11) 


-0.019(5) 


-0.003(8) 


0.018(5) 


0.017(8) 


-0.029(4) 


-0.013(6) 


0.025(3) 


0.021(6) 


-0.028(4) 


-0.014(4) 


0.024(3) 


0.018(4) 


-0.027(3) 


-0.014(4) 


0.023(3) 


0.017(3) 


-0.033(3) 


-0.019(3) 


0.027(2) 


0.019(3) 


-0.031(3) 


-0.019(3) 


0.026(2) 


0.018(3) 


-0.035(2) 


-0.023(2) 


0.029(2) 


0.020(2) 



mi 771,2 



0.02 





02 


0.02 





04 


0.02 





06 


0.02 





08 


0.04 





04 


0.04 





06 


0.04 





08 


0.06 





06 


0.06 





08 


0.08 





08 


0.02 





02 


0.02 





04 


0.02 





06 


0.02 





08 


0.04 





04 


0.04 





06 


0.04 





08 


0.06 





06 


0.06 





08 


0.08 





08 



-0.056(17) 

-0.065(14) 

-0.070(13) 

-0.074(12) 

-0.072(8) 

-0.075(7) 

-0.077(7) 

-0.082(6) 

-0.083(6) 

-0.090(5) 



(K + \(ud)s\p) 
0.067(48)0 
0.011(22)0 
0.014(16)0 
0.031(14)0 
0.038(11)0 
0.047(9) 
0.053(8) 
0.060(7) 
0.063(7) 
0.069(6) 



073(74) 

067(37) 

069(25) 

070(20) 

066(15) 

067(11) 

069(9) 

071(7) 

072(7) 

076(6) 



.068(19)0. 
.073(15)0. 
.076(13)0. 
.079(12)0. 
.081(8) 0. 
.082(7) 0. 
.083(7) 0. 
.088(6) 0. 
.088(6) 0. 
.093(5) 0. 
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TABLE IX: Continued from Table IVTTT1 





RL operator 


LL operator 


mi ni2 


V — Pi 


P = P2 


P = Pi 


P — P2 






(K + \(ds)u\p) 




0.02 0.02 


n n^/i (q\ 
— u.uo^y j 


— 058('43"l - 


-0 051 (~\A\ 

VJ . VJU ± I ± L ± 1 


— U.UOD^04l ) 


0.02 0.04 




—0 038f99 N l - 


-0 054H 1 

VJ . \JO'-±\ ± ± 1 


— U.UOU^Z / J 


0.02 0.06 


— U.UO l\t) 


— U.UoO^lO ) - 


-U.Uo / {\J) 


— u.uou^ry j 


0.02 0.08 


— u.uoy ) 


— D 0^7(1 1) - 

vj . wo i i loj 


VJ . VJVJVJ I U j 


n dm f i ^ 


0.04 0.04 


—0.035(5) 


—0 033no"l - 


-0 056f6"l 

VJ . VJOVJ I VJ ) 


—0.043(11) 


0.04 0.06 


— U.UO 1 ^ J 


—0 034(8"> - 

VJ . VJOt: 1 (J 1 


-0 058f5 N l 

VJ . UUUl (J J 


— U.U4: ( J 


0.04 0.08 




—0 036(7) - 

VJ • VJO VJ 1 1 1 


-0 060(5) 

VJ • VJVJVJ 1 U 1 


— U.UOU^ / J 


0.06 0.06 


— U.UoO^O j 


— 033f5") - 

U.UOOl O 1 


-0 061 f4"l 

vj . vjvj ± \ J 


— U.UDU^O J 


0.06 0.08 


U.UOOl O J 


—0 036(5) - 

\J • UOU 1 U ) 


-0 062(4) 

VJ . VJVJ^j V J 


— 053(5"! 


0.08 0.08 


— U.UO l j 


— 035(4) - 


-0 065(4) 

VJ . VJ U O \ t: 1 














0.02 0.02 


U.UlO^ / j 


079f49 N l 


069H 7") 

VJ . VJUii 11 1 1 


u.u / y ^oo j 


0.02 0.04 


n ni i ( K\ 


U.UOl) I ± cM 


U.UUvJl lO 1 


U.UOl^oZ J 


0.02 0.06 


n ni r\( K\ 

U.UIU^O ) 


0.023(13) 


0.070(11) 


U.UDl ^Zo J 


0.02 0.08 


u.uuy ) 


0.018(11) 


0.073(11) 


U.UDo^lo ) 


0.04 0.04 


010(3) 

V./ • V/ J. VJ l«J I 


0.016(8) 


0.068(7) 


053(13) 

VJ.VJ<J*Jl 1 'J i 


0.04 0.06 


0.010(3) 


0.014(6) 


0.071(6) 


0.057(10) 


0.04 0.08 


0.011(2) 


0.013(5) 


0.074(6) 


0.061(8) 


0.06 0.06 


0.009(2) 


0.010(4) 


0.075(5) 


0.061(7) 


0.06 0.08 


0.009(2) 


0.011(3) 


0.076(5) 


0.065(6) 


0.08 0.08 


0.008(1) 


0.009(2) 


0.079(4) 


0.067(5) 
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TABLE X: Fit results of the relevant form factor Wq with Eq. ()62|) for the principal matrix elements. 
Shown values are in lattice unit and by unrenormalized operators, dof = 16 for all. 





Co 


Cl 


C2 


C3 


x 2 


(7V°\(ud)RUL\p) 


-0.047(13) 


-0.125(31) 


-0.24(13) 


-0.27(8) 


14 (11) 


{ir°\(ud) L u L \p) 


0.067(15) 


0.070(33) 


0.09(16) 


0.09(10) 


0.6 (22) 


(K°\(u S ) R u L \p) 


0.055(12) 


0.071(32) 


0.07(14) 


0.12(7) 


1.6 (24) 


(K°\(u S ) L u L \p) 


0.026(9) 


0.023(23) 


-0.07(10) 


0.19(6) 


4.5 (61) 


(K+\(us) R d L \p) 


-0.021(7) 


-0.057(16) 


-0.13(7) 


-0.02(5) 


4.7 (55) 


(K + \(us) L d L \p) 


0.023(6) 


0.030(14) 


0.07(7) 


-0.02(4) 


1.9 (39) 


(K+\(ud) RSL \p) 


-0.046(14) 


-0.116(34) 


-0.18(15) 


-0.35(9) 


15 (12) 


(K + \(ud) L s L \p) 


0.072(15) 


0.069(35) 


0.06(16) 


0.15(11) 


0.4 (19) 


(K+\(ds) R u L \p) 


-0.033(8) 


-0.010(22) 


0.06(9) 


-0.10(5) 


0.5 (16) 


(K+\(ds) L u L \p) 


-0.049(11) 


-0.043(26) 


0.01(12) 


-0.17(7) 


0.5 (10) 


(r)\(ud) R u L \p) 


0.013(6) 


-0.012(15) 


-0.04(6) 


-0.02(4) 


4.7 (66) 


(i]\(ud) L u L \p) 


0.061(13) 


0.053(32) 


-0.01(14) 


0.20(9) 


0.5 (10) 



5(3 



